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A  method  is  proposed  which  aims  at  enhancing  the  performance  of  gen¬ 
eral  classes  of  elements  in  problems  involving  strain  localization.  The  method 
exploits  information  concerning  the  process  of  localization  which  is  readily 
available  at  the  element  level.  A  bifurcation  analysis  is  used  to  determine  the 
geometry  of  the  localized  deformation  modes.  When  the  onset  of  localization  is 
detected,  suitably  defined  shape  functions  are  added  to  the  element  interpola¬ 
tion  which  closely  reproduce  the  localized  modes.  The  extra  degrees  of  freedom 
representing  the  amplitudes  of  these  modes  are  eliminated  by  static  conden¬ 
sation.  The  proposed  methodology  can  be  applied  to  2D  and  3D  problems 
involving  arbitrary  rate-independent  material  behavior.  Numerical  examples 
demonstrate  the  ability  of  the  method  to  resolve  the  geometry  of  localized 
failure  modes  to  the  highest  resolution  allowed  by  the  mesh. 
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A  FINITE  ELEMENT  METHOD  FOR  LOCALIZED  FAILURE  ANALYSIS 


1.  Introduction 

A  phenomenon  that  frequently  accompanies  inelastic  deformation  is  the  formation  of  lo¬ 
calized  bands  of  intense  straining.  This  phenomenon  of  localization  occurs  in  a  wide  variety  of 
solids;  in  ductile  single  crystals.  Chang  and  Asaro  [l] ,  in  structural  metals,  Costin  et  al.  [2], 
Cox  and  Low  [3],  in  saturated  clays,  Vardoulakis  [4,5],  in  rocks,  Waversik  and  Brace  [6]  and 
in  concrete,  van  Mier  [7].  The  physical  mechanisms  responsible  for  triggering  localization  vary 
widely.  For  example,  at  high  rates  of  loading  in  structural  metals,  thermal  softening  plays  a 
key  role  in  initiating  localization  [2].  Localization  also  occurs  in  metals  at  low  rates  of  load¬ 
ing.  where  thermal  effects  are  negligible.  [1,3; .  Once  localization  takes  place,  large  strains  can 
accumulate  inside  the  band  and  lead  to  fracture. 

It  is  important  to  note  that  a  concentration  of  deformation  into  a  more  or  less  well  defined 
band  does  not  necessarily  constitute  localization  in  the  sense  used  here.  We  use  the  term 
localization  to  refer  to  situations  in  which  the  concentration  of  deformation  into  a  band  emerges 
as  an  outcome  of  the  constitutive  behavior  of  the  material.  Accordingly,  the  orientation  of  the 
band  is  characteristic  of  the  material,  rather  than  a  consequece  of  boundary  conditions.  We  will 
refer  to  such  a  band  of  localized  deformation  as  a  shear  band,  although  in  general  dilatation  as 
well  as  shear  takes  place  within  the  band. 

Since  large  strains  accumulate  inside  a  shear  band  without  substantially  affecting  the  strains 
in  the  surrounding  material,  e.  g.  Hutchinson  and  Tvergaard  ‘8  and  Saje  et  al.  [9],  shear 

bands  can  only  be  accurately  described  by  a  conventional  displacement  finite  element  method 

• 

if  the  band  interfaces  follow  element  boundaries.  Furthermore,  the  mesh  sets  the  minimum 
band  thickness  at  one  element  width.  Thus,  in  order  to  accurately  resolve  a  narrow  shear 
band,  small  elements  are  required,  with  element  boundaries  that  follow  band  directions.  Finite 
element  analyses  of  shear  bands  based  on  fine  meshes  of  quadrilateral  elements  built  up  from  four 
constant  strain  triangles  have  given  sharply  localized  deformation  modes,  see  e.  g.  Tvergaard 
et  al.  [10  ,  Tvergaard  [11[  and  the  review  of  Needleman  and  Tvergaard  12,.  A  material 
instability  analysis.  Hadamard  [13].  Hill  il  l  .  Mandel  [15,.  Thomas  16'.  Rice  17  ,  is  used  to 
orient  the  quadrilaterals.  Localization  is  found  in  numerical  calculations  even  when  the  mesh 


in  not  optimally  designed,  as  discussed  in  [10],  but  with  a  certain  delay  and  with  some  mesh 
induced  shear  band  broadening. 

In  a  two  dimensional  rectangular  mesh,  quadrilaterals  built  up  of  four  crossed  triangles  can 
resolve  narrow  bands  of  deformation  in  four  directions;  parallel  with  the  sides  of  the  quadrilat¬ 
eral  and  with  the  element  diagonals.  By  way  of  contrast,  isoparametric  quadrilaterals  can  only 
resolve  such  narrow  bands  parallel  with  the  element  sides.  Finite  element  analyses  of  oblique 
localized  deformation  bands  using  isoparametric  quadrilateral  elements,  e.  g.  Prevost  and 
Hughes  [18],  Prevost  [19],  Wiliam  et  al.  [20],  do  exhibit  the  tendency  for  shear  bands  to  form 
when  appropriate  critical  conditions  are  reached  but  have  not  exhibited  the  sharp  localization 
obtained  using  crossed  triangles. 

However,  the  crossed  triangle  finite  element  formulation  is  specifically  geared  to  two  di¬ 
mensional  problems  and  it  does  require  careful  mesh  design  to  account  for  likely  directions  of 
localization.  In  this  paper  a  method  is  proposed  which  aims  at  enhancing  the  performance 
of  isoparametric  elements  in  problems  involving  strain  localization.  The  formulation  here  is 
restricted  to  small  displacement  gradient  theory  and  to  rate  independent  constitutive  relations. 
The  method  uses  a  material  instability  bifurcation  analysis  [13-17],  carried  out  at  the  element 
level,  to  determine  the  geometry  of  the  ensuing  localized  deformation  modes.  When  the  onset 
of  localization  is  detected,  the  element  interpolation  is  extended  by  adding  to  it  suitably  de¬ 
fined  shape  functions  that  reproduce  the  localized  deformation  modes.  These  additional  shape 
functions  render  the  element  incompatible,  but  we  show  that  the  element  satisfies  the  patch 
test,  Irons  [21],  Taylor  et  a!.  [22],  The  extra  degrees  of  freedom  representing  the  amplitudes 
of  the  localized  modes  are  eliminated  by  static  condensation.  It  should  be  emphasized  that 
all  these  operations  can  be  implemented  at  the  element  level.  The  proposed  methodology  can 
be  applied  to  two  and  three  dimensional  problems  and  does  not  require  a  mesh  design  to  take 
explicit  account  of  likely  directions  of  localization.  Pietruszcak  and  Mroz  23  previously  incor¬ 
porated  a  shear  band  mode  of  deformation  into  a  finite  element  formulation.  However,  in  their 
approach,  the  localized  shear  band  is  specified  by  the  consitutive  relation  as  the  sole  mode  of 
inelastic  deformation  and  the  element  strain  rates  are  directly  given  in  terms  of  conventional 
shape  functions 
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Numerical  results  are  presented  in  Section  4  for  localization  in  both  nearly  incompressible 
solids  and  in  compressible  solids.  Even  in  the  absence  of  shear  bands,  isoparametric  elements 
are  of  limited  value  for  nealy  incompressible  solids,  such  as  the  classical  plastic  solid.  Mesh 
locking  occurs  as  a  result  of  the  near  incompressibilit  y  in  the  fully  plastic  range,  Nagtegaal  et  al. 
[24].  This  difficulty  can  be  overcome  through  the  use  of  a  method  to  handle  interna)  constrains; 
for  example,  here  we  use  a  B  technique,  Hughes  [25].  It  bears  emphasis,  however,  that  to  resolve 
localized  deformation  patterns  in  a  nearly  incompressible  solid  it  is  necessary,  but  not  sufficient, 
to  prevent  mesh  locking.  Our  numerical  results  show  that  the  enhanced  element  introduced 
here  does  resolve  localization  more  sharply  than  the  underlying  B-isoparametric  element.  Also, 
for  compressible  solids,  where  mesh  locking  due  to  incompressibility  is  not  an  issue,  numerical 
results  are  presented  that  illustrate  the  improved  performance  of  the  enhanced  element. 

2.  Localization  of  as  a  Bifurcation  Phenomenon 

In  this  section  we  review  some  aspects  of  the  general  theory  of  localization  of  inelastic 
deformations.  Some  of  the  basic  principles  underlying  the  theory  follow  from  Hadamard's 
studies  of  elastic  stability  [13],  extended  to  the  inelastic  context  by  Thomas  [16],  Hill  [14], 
Mandel  [15  and  Rice  [17] .  Here,  for  simplicity  attention  is  confined  to  infinitesimal  deformations 
and  thermally  decoupled,  rate-independent  materia!  behavior. 

Consider  a  homogeneous,  homogeneously  deformed  solid  subjected  to  quasi-static  incre¬ 
ments  of  deformation  <.  We  wish  to  determine  if  a  bifurcation  can  occur  in  such  a  manner  that 
subsequent  deformations  become  discontinuous  across  a  plane  of  orientation  n.  Let  u  be  the 
displacement  field  in  the  solid.  W'hereas  u  itself  remains  continuous  after  the  onset  of  local¬ 
ization,  the  displacement  gradients  Vu  will  exhibit  a  jump  across  the  plane  of  discontinuity,  i. 
e., 

K.jI  =  i=  0  (1) 

where  the  superindex  V!  refers  to  the  plus  side  of  the  plane  of  discontinuity  and  '  to  the 
minus  side.  Maxwell  s  compatibility  conditions  necessitate  that  the  jump  (1)  be  of  the  form 

=  (2) 
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for  some  vector  g.  Let  us  further  define  m  to  be  the  unit  vector  along  g,  i.  e., 

™ i  =  9,/ 9,  9=\  g!  (3) 

The  pair  of  unit  vectors  m  and  n  entirely  define  the  nature  of  the  discontinuity.  The  corre¬ 
sponding  strain  jump  is  given  by 

fc.i  ]  =  (s.  nj  +  9j  n,  )/2  (4) 

Often,  two  parallel  planes  of  discontinuity  pair  up  to  form  a  band  of  instensely  deformed 
material.  Two  limiting  cases  are  noteworthy: 

(a)  m  orthogonal  to  n.  The  material  in  the  band  deforms  in  simple  shear,  i.  e,  a  shear  band 
develops. 

(b)  m  parallel  to  n.  The  band  undergoes  extension  normal  to  the  planes  of  discontinuity.  In 
certain  circumstances  this  can  be  interpreted  as  a  splitting  failure  mode. 

Shear  bands  are  characteristic  of  materials  exhibiting  isochoric  plasticity.  On  the  other  hand, 
splitting  failure  modes  can  be  thought  of  as  idealizations  of  the  separation  processes  which  occur 
during  progressive  brittle  failure.  In  between  these  two  extremes  lies  a  continuous  spectrum  of 
mixed  failure  modes  for  which  m  and  n  are  neither  orthogonal  nor  parallel.  Plastic-fracturing 
materials  such  as  concrete,  rocks  and  ceramics  are  capable  of  failing  in  mixed  modes,  with  the 
angle  between  m  and  n  lying  anywhere  from  C f  to  9(7  [17,26.27!. 

We  next  investigate  under  what  conditions  localized  failure  modes  are  possible.  To  this 
end,  let  us  assume  that  the  solid  is  at  the  onset  of  localization.  At  this  point,  the  stresses  o 
and  strains  c  are  continuous  throughout  the  body.  However,  the  stress  and  strain  rates  a  and 
e  will  exhibit  discontinuities  across  a  plane  whose  orientation  is  to  be  determined.  Assuming 
rate-independent  material  behavior  the  incremental  stress-strain  relations  take  the  form 

—  Dijklt-kl  (5) 

where  D  is  the  tangent  stiffness  tensor  for  the  material.  Rate  independence  implies  that  D  is 
homogeneous  of  degree  zero  inc.  We  restrict  attention  to  cases  where  D  is  piecewise  independent 


of  c.  For  example,  in  classical  plasticity  D  has  two  branches,  one  corresponding  to  plastic 
loading,  the  other  to  elastic  unloading. 

For  elastic-plastic  solids,  the  localization  bifurcation  analysis  is  carried  out  using  the  same 
branch  of  the  moduli  D  across  the  incipient  planes  of  discontinuity,  which  corresponds  to  inves¬ 
tigating  bifurcation  for  a  linear  comparison  solid,  Hill  [33j,  Raniecki  and  Bruhns  [34].  Hence, 
taking  jumps  in  (5)  leads  to 

faj  j  =  Ay«  ¥ki  ]  (6) 

Equilibrium  across  the  discontinuity  planes  requires  that  the  tractions  t  be  continuous,  i.  e., 

IU  =  =  n,  [jdtJ|  =  °  (7) 

Combining  (C)  and  (7).  it  follows  that 

ft.  D.jki  f iki  {  =  0  (8) 

Finally,  using  the  kinematic  relation  (4)  and  the  definition  (3)  we  obtain 

Ajk{n)mk  =  (n,DtJk,m)mk  =  0  (9) 

This  condition  has  to  be  satisfied  by  m  and  n  for  the  localized  mode  to  be  possible.  The  onset 
of  localization  occurs  at  the  first  point  in  the  deformation  history  for  which  a  nontrivial  solution 
of  (9)  exists.  For  localization  to  occur  along  the  direction  n,  the  localization  matrix  A(n)  has 
to  have  at  least  one  zero  eigenvalue.  This  in  turn  necessitates 

/( n)  =  det(A(uj)  =  0  (10) 

If  a  unit  vector  n  satisfying  (10)  can  be  found,  the  corresponding  vector  m  which  completes  the 
definition  of  the  localized  mode  follows  from  (9). 

A  numerical  procedure  for  computing  m  and  n  from  the  above  equations  is  given  in  Ap¬ 
pendix  1.  In  two  dimensions,  the  localization  angles  can  be  computed  directly  as  the  roots  of  a 
cubic  polynomial,  while  in  three  dimensions  the  localization  directions  are  computed  iteratively. 
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An  example  of  application  of  the  solution  procedure  to  complex  three-dimensional  constitutive 
equations  is  found  in  [27]. 

3.  An  Assumed  Strain  Finite  Element  Method  for  Localized  Failure  Analysis 

It  has  been  long  recognized  [10],  that  isoparametric  elements  face  inherent  limitations  in 
the  presence  of  localized  deformations.  The  essence  of  the  problem  can  be  illustrated  by  means 
of  Fig.  1,  which  shows  a  shear  band  traversing  a  regular  mesh  at  45° .  The  elements  on  the 
boundary  of  the  shear  band  should  exhibit  a  strain  discontinuity.  However,  such  a  deformation 
pattern  is  not  well  represented  within  the  isoparametric  interpolation.  On  the  contrary,  the 
element  tries  to  conform  to  the  deformation  field  by  averaging  the  deformation  on  both  sides  of 
the  discontinuity.  As  a  result,  the  true  discontinuous  behavior  may  be  smeared  out  over  several 
elements  or  in  some  cases  precluded  altogether. 

A  numerical  procedure  which  circumvents  these  difficulties  is  proposed  next.  For  simplicity, 
attention  is  confined  to  infinitesimal  deformations.  The  method  can  be  described  as  follows: 

(i)  For  every  element  in  the  mesh  a  localization  analysis  of  the  type  described  in  Section  2  is 
carried  out  at  the  reduced  quadrature  points,  where  the  stresses  and  other  state  variables  are 
computed  most  accurately  [35] .  For  instance,  in  4-node  and  9-node  isoparametric  elements 
localization  is  investigated  at  the  centroid  and  the  2x  2  Gauss  points,  respectively.  This 
analysis  is  repeated  at  every  step  of  the  solution  process  until  the  onset  of  localization  is 
detected  at  one  or  more  reduced  quadrature  points. 

(ii)  From  this  point  on,  suitable  deformation  modes  which  exactly  reproduce  the  discontinous 
deformation  patterns  are  added  to  the  element  or  elements  where  localization  is  detected. 
The  extra  degrees  of  freedom  are  eliminated  at  the  element  level  by  means  of  static  con¬ 
densation 

This  approach  takes  advantage  of  the  fact  that  considerable  information  concerning  the  local¬ 
ization  process  can  be  readily  obtained  at  the  element  level.  The  method  aims  at  utilizing  this 
information  to  enhance  the  performance  of  the  element  in  the  presence  of  localized  failure. 

3.1  Localized  deformation  modes 

We  envision  a  generic  element  in  which  the  displacement  field  is  interpolated  as 
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(11) 


where  u,„  are  the  nodal  displacements,  A'„  (x)  are  the  interpolation  functions  and  the  sum 
extends  to  all  AT  nodes  in  the  element.  A  typical  example  is  given  by  the  serendipity  famiiy  of  4 
to  9  node  isoparametric  elements.  Let  us  further  denote  by  {£ir,  r  =  1,  . . . ,  X  Q}  the  coordinates 
of  the  reduced  quadrature  points.  Here  the  subindex  r  ranges  over  all  NQ  reduced  quadrature 
points  in  the  element. 

Based  on  the  state  variables  computed  at  the  points  a  localization  analysis  of  the  type 
described  in  Section  2  can  be  carried  throughout  the  solution  process.  Let  us  assume  that  in 
the  generic  element  under  consideration  and  at  some  stage  of  the  solution  several  bifurcations 
have  been  detected  resulting  in  XL  localized  modes.  The  geometry  of  these  deformation  modes 
is  fully  determined  by  pairs  of  unit  vectors  { (m,.,  ,  n,„  ),o  =  1,  . . . .  XL)  representing  the  dis¬ 
placement  direction  and  the  normal  to  the  surface  of  discontinuity,  respectively.  Furthermore, 
let  £,  denote  the  reduced  quadrature  point  at  which  the  ath  localized  mode  has  been  detected. 

Localized  deformation  modes  exhibit  a  strain  discontinuity  along  the  surface  of  normal  rt, 
passing  through  £,  .  In  order  to  accommodate  these  deformation  modes,  we  extend  the  element 
interpolation  by  means  of  suitably  defined  shape  functions  which  closely  reproduce  the  localized 
deformation  patterns.  To  this  end.  for  each  localized  mode  we  define  the  functions 


K,  W  = 


/  m,„  [n.  •  (X-  t,  )j,  if  n,  •  (x-  i,  )  >  0; 


l  0, 


otherwise. 


(12) 


and 


K  W  = 


-  m,,,  in,  •  (x  -  i,  )|F  if  n,  ■  (X  -  i,  )  <  0; 


0,  otherwise.  v 

where  no  sum  is  implied  on  the  index  a.  These  functions  exhibit  a  jump  across  the  surface  of 
discontinuity 


UK.  )<j  t  =  UK  )u  I  =  nj" 


(14) 


Consequently,  any  convex  combination 


M.u  (x)  =  (1  -  A„  )M~  (x)  +  A„  AC  (x) 


(15) 


of  M ^  and  M~t  exhibits  the  same  jump  (14).  Thus,  (15)  defines  a  one-parameter  family  of 
possible  shape  functions  for  the  localized  modes.  Adding  these  shape  functions  to  the  element 
interpolation,  the  local  displacement  field  becomes 

N  NL 

%  (x)  =  Y2  Ar„  (x)  +  Y2  T.  (x)  (16) 

a  =  l  u  =1 

where  {*ju  ,  a  =  1,  ...,A'L}  are  the  amplitudes  of  the  localized  modes.  It  should  be  noted 
that  shape  functions  (15)  are  incompatible,  i.  e.,  in  general  do  not  satisfy  the  C"  continuity 
requirement  across  the  element  boundary  and  the  element  becomes  nonconforming. 

The  element  strain  field  is  computed  from  (16)  to  be 


S'  L 


(x)  =  (n^xJ-rUj.,  (x))  /2  =  («,„A,I.;(x)-tiJ„A'„.,  (x))  7„  ((A/,„  )tJ  (x)-r(M,„  )„  (x)) 


.1  =  1  .1=1 

(1") 

For  subsequent  derivations  it  proves  convenient  to  recast  this  expression  in  matrix  form.  Let  Uj 
denote  the  nodal  displacement  array  and  u;  the  collection  of  amplitudes  %,  of  the  additional 
incompatible  modes.  Then  (17)  can  be  expressed  in  the  more  compact  form 


e(x)  =  BiUi  +  B;U;  (18) 

It  follows  from  (12)  and  (15)  that  the  compatibility  matrix  IL  takes  the  form 


(Bb) 


«.7'» 


A„  (m,„  Tij„  -  m,„  n,„  )  2.  if  n,  •  (x  -  £,  )  >  0: 

-  (1  -  A„  )(m,„  n.j„  -  nij„  n,,,  ),'2,  if  it,  •  (x-  )  <  0. 


so  that  the  element  strain  field  exhibits  a  jump 


(19) 


(x)‘  =  7.,  nJt,  +  mj,,  n,„  )/2  (20) 

across  the  ath  surface  of  discontinuity.  Thus,  the  parameters  q„  represent  the  magnitude  of 
the  strain  discontinuity  associated  with  the  ath  localized  mode. 
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It  is  implied  in  (18)  that  in  the  absence  of  localized  modes  the  element  strain  field  takes 
the  form  «(x)  =  Bi  (x)^  ,  where  Bi  results  from  applying  the  strain-displacemen  t  relations  to 
the  conforming  interpolation  (11).  It  was  noted  in  (24’  that  in  applications  involving  isochoric 
perfect  plasticity  this  formulation  results  in  an  overly  stiff  response  due  to  the  poor  handling  by 
the  element  of  the  near-incompressibily  which  accompanies  extensive  plastic  flow.  A  number 
of  techniques  have  been  proposed  to  alleviate  the  locking  phenomena  associated  with  internal 
constraints  such  as  incompressibility.  In  this  work  we  adopt  Hughes’  B-method  [25  which  for 
the  4-node  quadrilateral  consists  of  replacing  the  compatibility  matrix  Bj  by  a  modified  B,  of 
the  form 


B,(x)h  B/' '  (x)  -  f  Bi"'  (x)  -  (1  -  e)B|"'  (x, )  (21) 

Here.  B/'r  and  B["'  are  defined  so  as  to  give  the  deviatoric  and  volumetric  parts  of  the  strain 
tensor  when  multiplied  by  Uj  ,  x,  is  the  centroid  of  the  element  and  c  is  a  small  stabilization 
parameter  which  is  included  to  suppress  spurious  pressure  modes  [28  .  For  the  9-node  element. 
B["'  is  computed  at  the  reduced  quadrature  points  and  extrapolated  to  the  rest  of  the  element 
by  means  of  suitable  shape  functions.  Hughes’  B-method  is  closely  related  to  selective  reduced 
integration  techniques  and  to  mixed  methods,  see  e.  g.  the  review  of  Hughes  et  al.  [29  . 

Substituting  Bi  in  (18)  by  the  modified  compatibility  matrix  B,  the  element  strain  field 
finallv  takes  the  form 


e(x)  =  B,  Ui  -  B;U; 


Let  us  define  the  stiffness  submatrices 


in  =  [  B7  DB,  dV 

Jil' 

hz  =  KJ,  =  /  B; 

J  u' 

i:?  =j  B7  DB -d 


DBdY 


K:?  =  j  B:DB:rfr 

where  Q  signifies  the  domain  of  the  element  and  D  are  the  tangent  stiffness  tensor  of  the 
material.  Since  the  localized  modes  represent  internal  degrees  of  freedom  of  the  element,  the 
corresponding  amplitudes  u  can  be  eliminated  by  static  condensation  30  to  obtain 
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Uj  =  -  K21  Kjj 


(24) 


and  the  effective  element  stiffness  becomes 

K  —  lfu  —  K12  (25) 

It  should  be  emphasized  that  all  the  matrix  operations  indicated  in  (25)  can  be  carried  out  at 
the  element  level. 

As  can  be  seen  from  (25),  in  the  present  method  the  stiffness  matrix  Kn  of  the  underlying 
conforming  element  is  modified  by  means  of  a  term  which  represents  the  net  effect  of  the  added 
localized  modes.  In  Section  4  it  is  shown  by  way  of  numerical  examples  that  the  enhanced 
element  is  more  conducive  to  localization  and  provides  sharper  results  than  those  obtained 
from  the  B-method  alone. 

3.2  Patch  test 

In  Section  3.1  it  was  noted  that  the  introduction  of  localized  deformation  modes  renders  the 
element  nonconforming.  To  insure  that  the  convergence  properties  of  the  underlying  conforming 
element  are  not  altered,  the  incompatible  modes  have  to  be  defined  in  such  a  way  that  the  patch 
test  is  satisfied.  This  test  determines  whether  or  not  the  element  consistently  reproduces  states 
of  constant  strain.  If  so,  convergence  takes  place  within  each  element,  e.  g.  Strang  and  Fix 
128]. 

The  definition  of  shape  functions  (15)  contains  a  set  of  free,  as  yet  undetermined  parameters 
A„  .  Next,  it  is  shown  that  the  value  of  these  parameters  can  be  uniquely  determined  from  the 
requirement  that  the  element  satisfy  the  patch  test.  Taylor  et  al.  [22]  derived  a  necessary 
condition  for  the  patch  test  to  be  passed  by  a  nonconforming  element.  With  the  present 
notation  such  condition  reads 


In  view  of  (19),  (26)  necessitates 


BjdV  =  0 


(26) 


(-  (1  -  A„  )A'"  +  A„  A"*  j(m,„  tij 


+  "ij..  )/2  =  0 


(27) 


II 


where  A '+  and  A',  are  the  areas  of  the  element  subdomains  lying  on  the  plus  and  the  minus 
sides  of  the  ath  surface  of  discontinuity,  respectively.  Condition  (27)  is  satisfied  by  setting 


which  yields 


-(i  -  a„  )a-j  t\u  a-;;  =  0 


(28) 


(29) 


where  Ar  —  A,*  +  A't~  is  the  total  area  of  the  element.  In  practice,  satisfactory  results  are 
obtained  by  computing  the  areas  A'*  and  A~~  by  numerical  integration. 

3.3  Formulation  as  an  assumed  strain  method 

A  broad  class  of  finite  element  schemes  commonly  referred  to  as  6-methods  consist  of 
postulating  a  relation  <  =  Bu  for  some  suitably  defined  operator  B  and  a  tangent  stiffness 
matrix  of  the  form 


K  =  /  Br  D&dV  (30) 

J  ir 

This  class  of  methods  originated  with  the  pioneering  paper  of  Hughes  i25j  and  have  been  par¬ 
ticularly  successful  in  overcoming  mesh  locking  difficulties  associated  with  internal  constraints. 
The  formulation  proposed  here  can  be  rephrased  as  a  B-method  by  noting  that  (25)  can  be 
alternatively  expressed  as 

K  =  j  (6,  -  K13  B,  )r  D(B,  -  K13  K,1  B;  )dV  (31) 

J 

Comparison  of  (30)  and  (31)  warrants  the  identification 

B=  Bj  -  K1:  KJ;1  B.  (32) 

Thus,  the  proposed  formulation  can  be  regarded  as  a  B-method  in  which  the  the  compatibility- 
operator  is  postulated  to  take  the  form  (32). 

The  expression  (31)  for  the  effective  stiffness  matrix  can  be  directly  obtained  from  the 


Hellinger-Reissner  variational  principle.  This  principle  characterizes  the  equilibrium  and  com¬ 
patibility  equations  as  the  Euler  equations  of  the  potential 


L(u,c)  =  J  | Ajfci«. >(«.. J  +  -  Aj/c/f.jt/ci  -  /•u.J  rfV  -  J  t,u,ds  (33) 

where  f)  is  the  domain  occupied  by  the  body,  3  is  the  traction  boundary  and  f  and  t  denote  the 
body  forces  and  applied  tractions,  respectively.  Guided  by  the  results  obtained  from  the  static 
condensation  procedure  we  postulate  the  following  two-field  interpolation  for  displacements  and 


strains 


«.■(*)  =!>-*«(*)  (34) 

«(x)  =  Buj  =  (Bi  -  K12  K;2j  B2  )nL 

where,  as  before,  u,  denotes  the  collection  of  nodal  displacements  u,„  and  A'„  are  the  conformal 
shape  functions.  Substituting  the  interpolated  fields  (34)  into  the  Hellinger-Reissner  potential 
the  discretized  Euler  equations  are  found  to  be 


E  [£  BT  DBrfv]  u,  =  £  f, 

[  BrD(Bi  -  B)dV  Uj  =  0 

Jw 


where  the  element  forces  are  given  by 


(5  ).„  =  /  IN„dV  +  l  t,  A’„  dS 
Jur  -'an:' 


It  is  interesting  to  note  that  only  the  conformal  shape  functions  are  involved  in  the  computation 
of  the  element  forces.  It  follows  from  (35a)  that  the  element  stiffness  matrix  takes  the  form 


=  /  BT 

J  n« 


DBdE 


On  the  other  hand,  satisfaction  of  eq.  (35b)  for  arbitrary  incremental  displacements  necessitates 


f  BrD(B1  -  B)dV  =  0 
J  u- 


If  this  condition  is  identically  satisfied  the  stiffness  matrix  (37)  can  be  alternatively  expressed 
in  B-form 


K  =  /  Br  DBd  V  (39) 

J  w 

A  simple  computation  reveals  that  for  the  present  choice  of  B-matrix  (35b)  condition  (38)  is  in¬ 
deed  identically  satisfied.  Hence  the  method  can  be  directly  derived  from  the  Hellinger-Reissner 
variational  principle.  The  conditions  under  which  a  B-method  is  variationally  consistent  were 
first  elucidated  in  [31]  where  orthogonality  conditions  (38)  were  derived  within  the  more  general 
context  of  the  Hu-Washizu  principle. 

Satisfaction  of  orthogonality  conditions  (38)  not  only  places  the  method  on  a  solid  varia¬ 
tional  foundation  but  also  has  a  consequence  of  some  practical  importance.  Eq.  (35a)  indicates 
that  the  internal  forces  evolve  at  a  rate 


f7f  =  I  Bf DB</r 


Integrating  in  time  we  obtain  the  expression 


r  =E/r  B>dv 


Thus,  the  internal  forces  can  be  computed  by  using  the  conformal  compatibility  matrix  B; 
regardless  of  whether  localized  deformation  modes  have  been  added  to  the  interpolation.  On 
one  hand,  this  has  the  effect  of  eliminating  the  overhead  involved  in  computing  the  enhanced 
B-matrix  during  the  computation  of  the  internal  forces.  Secondly,  eq.  (41)  guarantees  that 
as  long  as  the  stress  history  is  continuous  in  time  the  internal  forces  will  also  be  a  continuous 
function  of  time.  In  particular,  no  sudden  jumps  will  be  experienced  by  the  internal  forces  on 
adding  a  localized  mode  to  the  interpolation. 

4.  Numerical  Examples 

Considerable  insight  into  the  performance  of  the  method  can  be  obtained  by  examining 
the  behavior  of  a  single  element.  A  first  simple  example  of  this  nature  is  shown  in  Fig.  2. 
The  problem  concerns  one  element  in  which  three  nodes  are  constrained  in  both  directions 
while  displacements  are  prescribed  on  the  fourth  node  which  result  in  nonuniform  shear  across 
the  element.  The  computed  force-displacement  diagrams  corresponding  to  an  isoparametric 
element,  Hughes'  B-method  and  the  present  approach  are  shown  in  Fig.  2.  It  is  seen  that  the 


*  /•  /•  /•«**\,*V*  .  V  *  •’*  .\V  /• 
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isoparametric  element  exhibits  considerable  residual  stiffness  once  the  element  is  fully  plastified. 
By  contrast,  both  the  B  method  and  the  enhanced  element  predict  limit  loads  of  similar  values. 
Once  these  limit  loads  are  attained  no  residual  stiffness  is  observed  and  a  failure  mechanism  is 
formed.  However,  the  nature  of  the  failure  mechanisms  predicted  by  both  methods  is  signifi¬ 
cantly  different.  Fig.  4  shows  the  computed  strain  histories  at  all  quadrature  points.  It  is  seen 
that  those  obtained  from  the  6-method  show  essentially  proportional  growth  and  no  tendency 
towards  localization.  By  contrast,  in  the  enhanced  element  a  bifurcation  is  detected  as  the 
centroid  plastifies  which  results  in  the  localized  mode  shown  in  Fig.  3.  From  this  point  on,  the 
strain  within  the  lower  half  of  the  element  remains  constant  and  all  subsequent  deformations 
are  localized  to  the  upper  half.  Thus,  whereas  adequate  handling  of  incompressibility  suffices 
to  capture  limit  loads  accurately  it  appears  that  to  obtain  sharp  predictions  of  localized  failure 
further  structure  has  to  be  added  to  the  element. 

A  second  example  is  concerned  with  shear  band  formation  in  a  perfectly  plastic  von  Mises 
solid  with  a  periodic  array  of  inhomogeneities,  Fig.  5.  The  body  is  subjected  to  an  overall 
uniform  shear  deformation.  The  inhomogeneities  are  modelled  by  one  element  voids  which  act 
as  nucleation  sites  from  where  shear  bands  emanate.  The  material  parameters  used  in  the 
analysis  are  E  =  2  x  10"' ,  v  —  0.3  and  av  -  200,  where  E ,  v  and  ay  are  the  Young's  modulus, 
Poisson’s  ratio  and  yield  stress,  respectively.  Fig.  6  shows  the  deformed  mesh  corresponding  to  a 
prescribed  overall  shear  deformation  7  =  1.6x  10~3  .  A  shear  band  spanning  the  inhomogeneities 
is  clearly  apparent.  The  bifurcated  elements  and  localization  directions  are  shown  in  Fig.  7. 
It  is  seen  that  localization  is  confined  to  a  diagonal  band  two  elements  wide.  Fig.  8  shows 
the  distribution  of  equivalent  plastic  strain  across  the  center  of  the  shear  band  computed  from 
the  B-method  and  the  enhanced  element  at  two  stages  of  the  solution  process.  Again  it  is 
seen  that  the  enhanced  element  results  in  a  sharper  shear  band  with  higher  amounts  of  plastic 
deformation  concentrated  over  a  narrower  region  than  the  B-method. 

Our  next  example  concerns  a  plane  strain  rigid  punch  problem  and  aims  at  assessing  the 
performance  of  the  method  in  the  presence  of  curved  shear  bands.  Fig.  9  shows  the  geometry 
of  the  problem  and  the  final  deformed  mesh  for  a  prescribed  displacement  \i  —  0.4.  A  perfectly 
plastic,  von  Mises  model  was  assumed  with  material  constants  E  -  30000,  v  =  0.3  and  cy  —  60. 


Fig.  10  shows  the  contours  of  maximum  shear  strains  computed  from  the  B-method  and  the 
enhanced  element.  A  well  defined  shear  band  is  observed  in  both  cases.  However,  the  computed 
shear  strains  are  seen  to  be  about  higher  in  the  case  of  the  enhanced  element.  Thus,  the 
proposed  method  is  seen  to  improve  the  performance  of  the  underlying  element  regardless  of 
the  orientation  of  the  shear  band  with  respect  to  the  mesh. 

The  examples  discussed  above  show  a  purely  quantitative  difference  in  the  response  as 
computed  with  or  without  localized  modes.  Our  final  example  is  concerned  with  a  case  in 
which  the  enhanced  element  accurately  captures  a  failure  mode  which  the  underlying  element 
fails  to  detect  Here  we  consider  a  rectangular  plane  strain  specimen  compressed  between  two 
rigid  platens.  Fig  11  Perfect  stick  between  the  specimen  and  the  platens  is  assumed.  The 
material  model  used  in  the  calculations  is  Drucker-Prager  s  model  for  soils  ,32  .  with  a  friction 
angle  of  3ff  .  a  dilaianry  and<  of  l(t  and  a  cohesion  of  5.  Perfect  plasticity  is  assumed  and 
the  elastic  response  is  takei.  to  b<  linear  and  isotropic  with  E  —  5000  and  r/  =  0.3.  This 
constitutive  relation  falls  within  the  framework  analyzed  by  Rudnicki  and  Rice  i26  ,  and  a 
noteworthy  outcome  of  their  instability  analysis  is  that  the  n unassociated  character  of  the 
constitutive  response  has  >•  «•!!*-«  t  of  promoting  localization  It  should  also  be  emphasized 
that  the  material  under  cun  i.n  t  'ioi,  i-  plastic  all\  dilatant  and  mesh  locking  due  to  near- 
lncompressibilit  y  is  of  n<  ton  rt;  I  ig  12  depico  the  load-displacement  curve  computed 
from  a  regular  isoparametric  element  and  the  same  element  with  localized  modes  added  to 
the  interpolation  As  can  be  seen  both  runes  are  virtually  identical  during  the  first  stages  of 
loading.  However,  ai  some  critical  deformation  the  enhanced  element  exhibits  a  sudden  drop 
in  the  load  This  corresponds  to  a  bifurcation  into  the  failure  mechanism  shown  in  Fig.  13b. 
By  contrast,  the  isoparametric  element  appears  to  preclude  the  development  of  a  failure  mode 
and  the  deformation  pattern  remains  essentially  unchanged  throughout  the  loading  process. 
Fig.  13a.  The  distribution  of  effective  plastic  strain  shown  in  Fig  14  further  illustrates 
this  discrepancy.  Whereas  the  isoparametric  element  predicts  a  smooth  distribution  of  plastic 
flow  over  the  specimen.  Fig.  14b.  the  enhanced  element  yields  a  rather  intricate  pattern  with 
localized  regions  of  intense  straining.  Fig  14a.  Since  in  this  example  incompressibilit  y  is  not  a 
issue,  it  is  clear  that  the  inability  of  the  isoparametric  element  to  produce  localized  failure  modes 


stems  from  reasons  other  than  mesh  locking.  As  discussed  in  Section  3,  the  smoothness  of  the 
isoparametric  interpolation  and  the  tendency  of  the  element  to  smear  out  strain  discontinuities 
are  at  the  basis  of  the  poor  performance. 

5.  Summary  and  Conclusions 

A  method  has  been  presented  which  aims  at  enhancing  the  performance  of  genera]  classes  of 
elements  in  problems  involving  localized  failure.  The  method  exploits  the  fact  that  considerable 
information  concerning  the  process  of  localization  can  be  readily  obtained  at  the  element  level. 
A  simple  analysis  suffices  to  detect  the  onset  of  localization  in  an  element  and  to  determine  the 
geometry  of  the  localized  deformation  modes.  This  information  is  utilized  to  set  up  additional 
shape  functions  which  closely  reproduce  the  localized  deformation  patterns.  The  additional 
degrees  of  freedom  are  then  eliminated  at  the  element  level  by  static  condensation.  Although 
the  resulting  element  is  nonconforming  the  patch  test  is  satisfied  and  the  convergence  properties 
of  the  underlying  conforming  element  are  not  altered.  The  method  can  be  formulated  as  a  B- 
method  or  derived  directly  from  the  Hellinger-Reissner  principle.  Numerical  examples  indicate 
that  whereas  limit  loads  can  be  accurately  computed  by  means  of  finite  element  methods  which 
adequately  handle  the  incompressibility  constraint,  to  obtain  sharp  predictions  of  localized 
failure  modes  using  quadrilateral  isoparametric  elements  further  structure  has  to  be  added 
to  the  interpolation.  The  present  approach  provides  an  efficient  means  of  incorporating  this 
additional  structure  into  general  classes  of  2D  and  3D  elements. 

Appendix  I.  Solution  of  the  Localization  Condition 

In  this  appendix  we  discuss  a  computational  procedure  for  detecting  the  onset  of  localization 
and  determining  the  localization  directions  m  and  a  It  proves  convenient  to  treat  the  2D  and 
2D  cases  separately . 

Three-dimensional  case.  If  the  material  response  is  path-dependent,  the  integration  of  the 
constitutive  equations  has  to  be  carried  out  incrementally.  If,  for  instance,  the  instantaneous 
initial  response  of  the  materia!  is  isotropic-elastic,  the  localization  function  (10)  initially  takes 
the  value 


/( n)  =  (A„  -2p,.)/r 


(.41.1) 


17 


where  A„  and  p„  are  the  Lame  constants  of  the  virgin  material.  Note  that  /( n)  >  0  and 
independent  of  n.  As  the  process  of  deformation  progresses,  /  may  develop  minima  which 
eventually  become  negative  or  zero.  This  in  turn  signals  the  onset  of  localization.  To  detect 
precisely  when  this  happens,  the  minima  of  the  localization  function  /  can  be  computed  at 
every  deformation  increment  and  the  localization  condition  (10)  checked  at  the  minima.  This 
leads  to  considering  the  constrained  minimization  problem 


minimize  /( n)  =  det(n,  D„ki  n,) 
subject  to  |  n  |  =  1 

where  D  is  the  current  value  of  the  tangent  moduli.  The  minima  are  characterized 
condition 


{A  1-2) 
by  the 


(A1.3) 


where  A  is  a  Lagrange  multiplier.  Differentiating  /  in  (A  1.3)  one  obtains 


Introducing  the  notation 


det(A.{n))D,  ,ki  Akjl  (n)n,  -  An,  =  0 


(A1.4) 


J'i( »)=  det(A(n))Dl]klAkjl  (n) 


the  minimum  condition  (A1.4)  can  be  recast  as 


Ml-5) 


J,i  (u)n,  -  An,  =  0  (A  1.6) 

The  solutions  of  (A  1.6)  can  be  found  in  two  steps. 

i)  Expressing  n  in  terms  of  spherical  angles,  i.  e.,  setting  n  =  (cosOcosS ,  cososinO ,  sind>). 
the  range  of  variation  0,  2x  \  x  [0,tt/2;  of  ( 9,o )  is  swept  at  5-degree  increments  to  determine  a 
first  approximation  n101  to  the  minima. 

ii)  The  locations  of  the  minima  are  then  pinpointed  by  means  of  the  iterative  scheme: 


r  /  I  k)  \  IV+11 


.  Ik-  I  )  Ik-  1  I 


Thus,  at  every  iteration  an  eigenvalue  problem  is  formulated  based  on  the  matrix  J  evaluated 
at  the  previous  iteration  n1*1 .  The  minimum  eigenvector  of  (A1.7)  is  taken  as  the  new  iterate 
n'*+1> . 

Once  the  orientations  of  the  discontinuity  planes  have  been  determined,  the  corresponding 
m- vectors  tire  computed  as  the  zero-eigenvectors  of  the  localization  matrix  A(n). 

Two-dimensional  case.  In  the  2D  case  the  localization  matrix  A(n)  is  2  x  2  and  one  readily 
finds  that 


det{  A(n))  =  d,  +  oj  n®  fa  +  a?  n*  +  dj  rij  t  ^ 

where 

a,,  =  Ain  A212  ~  A112A21) 
ai  =  A111A232  +  A111A212  —  A112A211 

On  —  A?222  112  ^1222  A  -^121 1  Al2 12 

(X3  —  ^1112^2222  -Dl21 1  Al222  Dj  122 -Al212 

dj  —  D 1212  —  D->2i2  D|222 

Setting  rq  =  cos0,  no  =  sinO  in  eq.  (A1.8)  the  localization  condition  becomes 

/(z)  =  a,,x4+ai3f5-!-aoz::-t-aoZ-!-a4=0  (41.10) 

where  one  writes  z  =  tanB .  In  general,  the  polynomial  /(x)  in  (A1.10)  is  positive  everywhere 
prior  to  localization.  Thus  the  onset  of  localization  can  be  determined  by  simply  examining 
the  sign  of  the  minima  of  /(x).  These  occur  at  the  roots  of  the  cubic  polynomial  f  (x)  which 
can  be  computed  in  close  form  by  means  of  Cardan’s  formulae.  As  long  as  the  minima  of  /(z) 
remain  positive  localization  does  not  develop.  The  onset  of  localization  is  signaled  by  one  of 
more  minima  of  /(x)  crossing  the  x-axis. 
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Figure  Captions 

Figure  1  Close  up  view  of  the  central  portion  of  the  shear  band  shown  in  Fig.  6. 

Figure  2  Single  element  test  problem. 

Figure  3  Force-displacement  curves  computed  from  an  isoparametric  element,  a  B  method  and 
the  enhanced  element. 

Figure  4  Evolution  of  shear  strains  in  test  problem  of  Fig.  2.  The  results  correspond  to  (a)  B 
method  and  (b)  enhanced  element. 

Figure  5  Unit  cell  for  a  periodic  distribution  of  voids  subjected  to  an  overall  shear  deformation. 
Figure  6  Deformed  mesh  for  the  problem  defined  in  Fig. 5  at  a  prescribed  shear  strain  7  =  0.0016. 
Magnification  factor  for  displacements  =  50. 

Figure  7  Distribution  of  bifurcated  elements  corresponding  to  the  deformation  shown  in  Fig.  6. 

The  arrows  represent  the  characteristic  directions  associated  with  the  localized  modes. 
Figure  8  Distribution  of  effective  plastic  strain  across  the  shear  band  at  two  stages  of  deforma¬ 
tion.  The  values  displayed  are  computed  at  the  element  centroids. 

Figure  9  Plane  strain  punch  problem  on  a  tapered  specimen.  Deformed  mesh  for  11  —  U.4. 
Figure  10  Contours  of  maximum  shear  strain  for  punch  problem  shown  in  Fig.  9  at  u  =  0.4 
Results  computed  from  (a)  B  method  and  (b)  enhanced  element. 

Figure  11  Plane  strain  soil  specimen  compressed  between  two  rigid  platens. 

Figure  12  Load-displacement  curve  for  problem  shown  in  Fig,  11.  Results  computed  from  (a) 
isoparametric  element  and  (b)  enhanced  element. 

Figure  13  Deformed  mesh  computed  from  (a)  isoparametric  element  and  (b)  enhanced  element. 

The  results  correspond  to  points  A  and  B  on  the  load-displacement,  diagram.  Fig.  12. 
Magnification  factor  for  displacements  =  25. 

Figure  14  Contours  of  effective  plastic  strain  computed  from  (a)  isoparametric  element  and  (b) 


enhanced  element. 
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